% This code can be executed to generate Figure S3. Raw data (load) is
% plotted for two samples: 25 and 15. A control tensile test of a single 
% tape is shown in green. The non-linearized elastic model (EF). is plotted
% as dotted-dashed lines. The linearized elastic model (EFLin) is plotted
% as dotted lines. The point where the elastic and peeling energies are 
% equal is marked with a yellow triangle (intersectUFelas). The point where 
% the elastic model becomes linear is marked with an orange square 
% (linearpts). The maximum loads are marked by grey dots (peakload).
% Additional formatting is done externally to exactly match the published 
% Figure S3.
% The folder titled, "Raw Data Test Files," must be open to run this code!
clear;
t=0.058; %mm
w=12.7; %mm
E=467.9745; % Pa
gamma=10.3014/1000;  %N/mm
beta=mean([189.393 190.444 184.519 182.807 186.996]); %mm
slope=(E*w)./1000; %N/mm

plotnorm=w;
%Determine which tensile data sets to plot
names=['2020-04-15_003_15.0.txt',...
    '2020-05-04_045_00.0.txt',...
    '2020-04-15_007_25.0.txt'];

dat_offset=[0.17 -0.52 0.14 0];
text_offset=[0.5 0.5 0];

%short algorithm to determine where each file name starts and ends to
%determine the junction angle information which is stored in the filenames
for i=1:length(names)
    chunk=i+3;
    test=names(i:chunk) == '.txt';
    if all(test)
        numchar=chunk;  % number of characters in each file name
        break
    end
end

%Initialize variables
ntrials=length(names)/numchar;
jj=zeros(1,ntrials);
jj(1)=1;
kk=zeros(1,ntrials);
kk(1)=numchar;
angles=zeros(1,ntrials);
peakloads=zeros(10,ntrials);
peakload=zeros(1,ntrials);
xmaxs=zeros(1,ntrials);
pklocs=zeros(1,ntrials);
peelmax=zeros(1,ntrials);
linearpts=zeros(1,ntrials);
intersectUx=zeros(1,ntrials);
intersectUFelas=zeros(1,ntrials);
Lo=zeros(1,ntrials);

colors=zeros(ntrials,2,3);
colors(1,1,:)=[70/255 0/255 130/255]; %data 1
colors(1,2,:)=[170/255 70/255 255/255]; %model 1
colors(2,1,:)=[0/255 120/255 80/255]; %data 2
colors(2,2,:)=[60/255 240/255 165/255]; %model 2
colors(3,1,:)=[70/255 150/255 180/255]; %data 3
colors(3,2,:)=[10/255 200/255 255/255]; %model 3

adjustloads=[];
adjustloc=[];
hold on;
for i=1:ntrials
    clear Data displacement load time
    filename=names(jj(i):kk(i));
    fileID=importdata(filename);
    j=jj(i)+numchar;
    k=kk(i)+numchar;
    jj(i+1)=j;
    kk(i+1)=k;
    Data=fileID.data;
    displacement=Data(:,1)-dat_offset(i);  %mm
    load=Data(:,2);  %N
    time=Data(:,3);  %s
    angles(i)=str2num(filename(16:19));
    Lo(i)=max(displacement);
    %stress and strain are not used, but they can be calculated as follows:
    %strain=(displacement./Lo(i)).*100;  %%
    %stress=((load)./Area); %MPa
    plot((displacement-dat_offset(i))./plotnorm,load,'LineWidth',2','Color',colors(i,1,:));
    
    phi=angles(i);
    xmax=@(phi) 2.*w.*tand(phi./2);
    fun11=@(x) 2.*gamma.*(x+dat_offset(4)).*((cscd(phi./2)).^2).*(cotd(phi./2));
    fun12=@(x) 2.*gamma.*(xmax(phi)-(x+dat_offset(4))).*((cscd(phi./2)).^2).*(cotd(phi./2));
    xxl1=linspace(0,xmax(phi)./2,1000);
    xxl2=linspace(xmax(phi)./2,xmax(phi),1000);
    y1=(fun11(xxl1));
    y2=(fun12(xxl2));
    
    %Optional: Plot Peeling Model
    %plot(xxl1./plotnorm,y1,'Color',colors(i,2,:),'LineWidth',2,'LineStyle','--');
    %plot(xxl2./plotnorm,y2,'Color',colors(i,2,:),'LineWidth',2,'LineStyle','--');
    text((xmax(phi)./2+text_offset(i))./plotnorm,((fun11(xmax(phi)./2)))+.2,[num2str(phi) '\circ'],'FontSize',14,'Color',colors(i,1,:));
    %Determine peaks of the data sets similarly to Figure 4A
    adjusttrials=[13 15 17 18 19 20 22 24]; 
    if any(ismember(adjusttrials,str2num(filename(12:14))))
        [pks, pkloc]=findpeaks(load,displacement);
        [pks2, pkloc2]=findpeaks(pks,pkloc);
        peaks=[pks2'; pkloc2']';
        sortpks=sortrows(peaks,1);
        peakloads(1:length(pks2),i)=sortpks(:,1);
        peakload(i)=peakloads(length(sortpks)-1,i);
        pklocs(i)=sortpks(length(sortpks)-1,2);
        adjustloads=[adjustloads max(pks)];
        adjustloc=[adjustloc sortpks(length(sortpks),2)-dat_offset(i)];
    else
        [pks, pkloc]=findpeaks(load,displacement);
        peaks=[pks'; pkloc']';
        sortpks=sortrows(peaks,1);
        peakload(i)=max(pks);
        pklocs(i)=sortpks(length(sortpks),2)-dat_offset(i);
    end
    xmaxs(i)=xmax(phi)./2;
    peelmax(i)=max(y1);
    
    %Determine & Plot Non-Linearized Elastic Model
    xx=linspace(0,xmax(phi),1000);
    EFfun=@(x) -(1/2).*E.*t.*x.*tand(phi).*log(abs((2.*x.*cotd(phi./2)-beta.*tand(phi))./(-beta.*tand(phi))));
    EF=EFfun(xx);
    plot(xx./plotnorm,EF,'Color',colors(i,2,:),'LineWidth',2,'LineStyle','-.');
    
    
    %Determine & Plot Linearized Elastic Model
    EFLin=EF;
    for n=1:(length(EFLin)-1)
        m=(EFLin(n+1)-EFLin(n))/(xx(n+1)-xx(n));
        b=EFLin(n)-m*xx(n);
        if m>slope
            xremain=linspace(xx(n),max(xx),(length(xx)-n+1));
            EFLin(n:length(xx))=xremain.*slope+b;
            break
        end
        linearpts(i)=xx(n);
    end
    plot(xx./plotnorm,EFLin,'Color',colors(i,2,:),'LineWidth',2,'LineStyle',':');

    %Determine Energy Intersect
    xmax=xmax(phi);
    Upeelinc=@(x) 2.*gamma.*((cscd(phi./2)).^2).*(1./2).*(x.^2).*(cotd(phi./2));
    Upeeldec=@(x) 2.*gamma.*((cscd(phi./2)).^2).*((1)./2).*((-1.*(xmax-x).^2)+2.*(xmax./2).^2).*(cotd(phi./2));
    xx=linspace(0,xmax,1000);
    xxl1=linspace(0,xmax./2,length(xx)./2);
    xxl2=linspace(xmax./2,xmax,length(xx)./2);
    xx=[xxl1 xxl2];
    Upeel=[Upeelinc(xxl1) Upeeldec(xxl2)];
    Uelas=cumtrapz(xx,EF);
    subU=Uelas-Upeel;
    for a=1:length(subU)-1
        test=subU(a);
        if test>0
            break
        end
    end
    intersectUx(i)=xx(a);
    intersectUFelas(i)=EF(a);

end
%Optional: Plot Peeling Maxima
%plot(xmaxs./plotnorm,peelmax,'Marker','.','Color',[180/255 0 0],'LineStyle','none','MarkerSize',30);

%Plot Max Load
scatter(pklocs./plotnorm,peakload,70,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0 0 0]...
    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',3/3);

%Plot Elastic Max 
scatter(adjustloc./plotnorm,adjustloads,70,'filled','MarkerFaceAlpha',3/3,'MarkerFaceColor',[255/255 139/255 130/255]...
    ,'MarkerEdgeColor',[180/255 68/255 54/255],'MarkerEdgeAlpha',3/3);

%Plot where Elastic Force meets alpha (orange square)
plot(linearpts./plotnorm,EFfun(linearpts),'Marker','s','MarkerFaceColor',[255/255 126/255 24/255],'LineStyle','none','MarkerSize',10);

%Plot Elastic Energy=Peeling Energy (yellow triange)
plot(intersectUx./plotnorm,intersectUFelas,'Marker','^','MarkerFaceColor',[255/255 193/255 37/255],'LineStyle','none','MarkerSize',10);

%Optional: Determine & plot Peeling Model
%xmaxfun=@(phi) w.*tand(phi./2);
%anglescontinuous=linspace(10,110,1000);
%peelmax=zeros(1,length(anglescontinuous));
%for j=1:length(anglescontinuous)
%    phi=anglescontinuous(j);
%    xmax=xmaxfun(phi);
%    fun11=@(x) 2.*gamma.*(x).*((cscd(phi./2)).^2).*(cotd(phi./2));
%    peelmax(j)=fun11(xmax);
%end
%plot(xmaxfun(anglescontinuous)./plotnorm,peelmax,'Color',[180/255 0 0],'LineWidth',2,'Marker','none');

%Formatting the plot
ylim([0 26]);
xlim([0 0.9])
xlabel('Tensile Displacement (mm)','FontSize',14);
ylabel('Load (N)','FontSize',14);
hold off;
ax=gca;
ax.FontSize=14;